home *** CD-ROM | disk | FTP | other *** search
/ Practical Algorithms for Image Analysis / Practical Algorithms for Image Analysis.iso / TARFILE.GZ / tarfile / ch_4.5 / xbdy / bdy.c next >
Encoding:
C/C++ Source or Header  |  1999-09-11  |  15.1 KB  |  602 lines
  1. /* 
  2.  * bdy.c
  3.  * 
  4.  * Practical Algorithms for Image Analysis
  5.  * 
  6.  * Copyright (c) 1997, 1998, 1999 MLMSoftwareGroup, LLC
  7.  */
  8.  
  9. /*
  10.  * BDY.C
  11.  *
  12.  * functions implementating C.T.Zahn's algorithm for boundary
  13.  * detection and representation
  14.  *
  15.  */
  16. #include <stdio.h>
  17. #include <malloc.h>
  18. #include <stdlib.h>
  19. #include <math.h>
  20. #include "ip.h"
  21. #include "bdy.h"
  22.  
  23. #define    WINDOW_ROWS    3
  24. #define NDIRECTIONS    8
  25.  
  26. #define    THRESH        1
  27. #define MIN_POLY_SIZE    1
  28. #define    ZERO        0
  29. #define ROOT2        1.414213562
  30.  
  31. #define    ON        1
  32. #define    OFF        0
  33. #undef    SHOW_CURV_PT
  34. #undef  DPRINT
  35.  
  36. #define    RESET        ON
  37.  
  38. #define IN_LIST        1
  39. #define    OUT_LIST    0
  40.  
  41.  
  42. /* global variables */
  43. extern struct curv_points *curv_head_in[NDIRECTIONS];
  44. extern struct curv_points *curv_head_out[NDIRECTIONS];
  45. extern struct curv_points *curv_tail_in[NDIRECTIONS];
  46. extern struct curv_points *curv_tail_out[NDIRECTIONS];
  47.  
  48. extern struct polygon *poly_head_ptr;
  49. extern struct polygon *poly_tail_ptr;
  50.  
  51. extern float *delta_phik, *delta_lk;
  52. extern long ncp;
  53.  
  54. unsigned int d1, d2, d3, d4, d5, d6, b1, b2, b3, b4, b5, b6;
  55. unsigned nbytes = 1;
  56.  
  57.  
  58. /*
  59.  * generate curvature points using Zahn algorithm
  60.  * (ref: C. T. Zahn, SLAC report No. 70, Dec. 1966)
  61.  */
  62. void
  63. curvature_points (unsigned char *window[], int ir, int jmax, int imax)
  64. {
  65.   int jc;
  66.  
  67.   for (jc = 0; jc < jmax - 2; jc++) {
  68.     d1 = *(*(window + (ir % WINDOW_ROWS)) + jc);
  69.     d2 = *(*(window + (ir % WINDOW_ROWS)) + jc + 1);
  70.     d3 = *(*(window + (ir % WINDOW_ROWS)) + jc + 2);
  71.     d4 = *(*(window + ((ir + 1) % WINDOW_ROWS)) + jc);
  72.     d5 = *(*(window + ((ir + 1) % WINDOW_ROWS)) + jc + 1);
  73.     d6 = *(*(window + ((ir + 1) % WINDOW_ROWS)) + jc + 2);
  74.     horizontal (ir, jc);
  75.   }
  76.   for (jc = 0; jc < jmax - 1; jc++) {
  77.     b1 = *(*(window + (ir % WINDOW_ROWS)) + jc + 1);
  78.     b2 = *(*(window + ((ir + 1) % WINDOW_ROWS)) + jc + 1);
  79.     b3 = *(*(window + ((ir + 2) % WINDOW_ROWS)) + jc + 1);
  80.     b4 = *(*(window + (ir % WINDOW_ROWS)) + jc);
  81.     b5 = *(*(window + ((ir + 1) % WINDOW_ROWS)) + jc);
  82.     b6 = *(*(window + ((ir + 2) % WINDOW_ROWS)) + jc);
  83.     vertical (ir, jc);
  84.   }
  85.   if (ir == imax - WINDOW_ROWS) {
  86.     ir++;
  87.     for (jc = 0; jc < jmax - 2; jc++) {
  88.       d1 = *(*(window + (ir % WINDOW_ROWS)) + jc);
  89.       d2 = *(*(window + (ir % WINDOW_ROWS)) + jc + 1);
  90.       d3 = *(*(window + (ir % WINDOW_ROWS)) + jc + 2);
  91.       d4 = *(*(window + ((ir + 1) % WINDOW_ROWS)) + jc);
  92.       d5 = *(*(window + ((ir + 1) % WINDOW_ROWS)) + jc + 1);
  93.       d6 = *(*(window + ((ir + 1) % WINDOW_ROWS)) + jc + 2);
  94.       horizontal (ir, jc);
  95.     }
  96.   }
  97. }
  98.  
  99.  
  100. /*
  101.  * check 3x2 horizontal array (in a binary image) and determine 
  102.  * direction codes for incoming and outgoing contour segment
  103.  */
  104. int
  105. horizontal (i, j)
  106.      int i, j;
  107. {
  108.   int in, out;
  109.   int err_flag = -1;
  110.  
  111.   if (d2 >= THRESH) {
  112.     if (d5 == ZERO) {
  113.       if (d6 >= THRESH)
  114.         in = 3;
  115.       else if (d3 >= THRESH)
  116.         in = 4;
  117.       else
  118.         in = 5;
  119.  
  120.       if (d4 >= THRESH)
  121.         out = 5;
  122.       else if (d1 >= THRESH)
  123.         out = 4;
  124.       else
  125.         out = 3;
  126.     }
  127.     else
  128.       return (err_flag);
  129.   }
  130.   else {
  131.     if (d5 >= THRESH) {
  132.       if (d1 >= THRESH)
  133.         in = 7;
  134.       else if (d4 >= THRESH)
  135.         in = 0;
  136.       else
  137.         in = 1;
  138.  
  139.       if (d3 >= THRESH)
  140.         out = 1;
  141.       else if (d6 >= THRESH)
  142.         out = 0;
  143.       else
  144.         out = 7;
  145.     }
  146.     else
  147.       return (err_flag);
  148.   }
  149.  
  150. /*
  151.  * curvature point found! 
  152.  * allocate memory and store it in linked list curv_points structure 
  153.  */
  154.   if (in != out) {
  155.  
  156. #ifdef SHOW_CURV_PT
  157.     printf ("   horizontal:  ir= %d, jc= %d, in= %d, out= %d\n",
  158.             2 * (i + 1), 2 * (j + 1) + 1, in, out);
  159. #endif
  160.     if (curv_head_in[in] == NULL) {
  161.       curv_head_in[in] = (struct curv_points *) calloc (nbytes, sizeof (struct curv_points));
  162.       curv_head_in[in]->prev = NULL;
  163.       curv_tail_in[in] = curv_head_in[in];
  164.     }
  165.     else {
  166.       curv_tail_in[in]->next = (struct curv_points *) calloc (nbytes, sizeof (struct curv_points));
  167.       curv_tail_in[in]->next->prev = curv_tail_in[in];
  168.       curv_tail_in[in] = curv_tail_in[in]->next;
  169.     }
  170.     if (curv_head_out[out] == NULL) {
  171.       curv_head_out[out] = (struct curv_points *) calloc (nbytes, sizeof (struct curv_points));
  172.       curv_head_out[out]->prev = NULL;
  173.       curv_tail_out[out] = curv_head_out[out];
  174.     }
  175.     else {
  176.       curv_tail_out[out]->next = (struct curv_points *) calloc (nbytes, sizeof (struct curv_points));
  177.       curv_tail_out[out]->next->prev = curv_tail_out[out];
  178.       curv_tail_out[out] = curv_tail_out[out]->next;
  179.     }
  180.     ncp++;
  181.  
  182.     curv_tail_in[in]->same_point = curv_tail_out[out];
  183.     curv_tail_in[in]->next = NULL;
  184.     curv_tail_out[out]->next = NULL;
  185.     curv_tail_in[in]->xn = curv_tail_out[out]->xn = 2 * (j + 1) + 1;
  186.     curv_tail_in[in]->yn = curv_tail_out[out]->yn = 2 * (i + 1);
  187.     curv_tail_in[in]->edge_in = curv_tail_out[out]->edge_in = in;
  188.     curv_tail_in[in]->edge_out = curv_tail_out[out]->edge_out = out;
  189.   }
  190.   return (1);
  191. }
  192.  
  193.  
  194. /*
  195.  * check 2x3 vertical array (in a binary image) and determine 
  196.  * direction codes for incoming and outgoing contour segment
  197.  */
  198. int
  199. vertical (i, j)
  200.      int i, j;
  201. {
  202.   int in, out;
  203.   int err_flag = -1;
  204.  
  205.   if (b2 >= THRESH) {
  206.     if (b5 == ZERO) {
  207.       if (b6 >= THRESH)
  208.         in = 1;
  209.       else if (b3 >= THRESH)
  210.         in = 2;
  211.       else
  212.         in = 3;
  213.  
  214.       if (b4 >= THRESH)
  215.         out = 3;
  216.       else if (b1 >= THRESH)
  217.         out = 2;
  218.       else
  219.         out = 1;
  220.     }
  221.     else
  222.       return (err_flag);
  223.   }
  224.   else {
  225.     if (b5 >= THRESH) {
  226.       if (b1 >= THRESH)
  227.         in = 5;
  228.       else if (b4 >= THRESH)
  229.         in = 6;
  230.       else
  231.         in = 7;
  232.  
  233.       if (b3 >= THRESH)
  234.         out = 7;
  235.       else if (b6 >= THRESH)
  236.         out = 6;
  237.       else
  238.         out = 5;
  239.     }
  240.     else
  241.       return (err_flag);
  242.   }
  243.  
  244. /*
  245.  * curvature point found! 
  246.  * allocate memory and store it in linked list curv_points structure 
  247.  */
  248.   if (in != out) {
  249.  
  250. #ifdef SHOW_CURV_PT
  251.     printf ("  vertical:  ir= %d, jc= %d, in= %d, out= %d\n",
  252.             2 * (i + 1) + 1, 2 * (j + 1), in, out);
  253. #endif
  254.  
  255.     if (curv_head_in[in] == NULL) {
  256.       curv_head_in[in] = (struct curv_points *) calloc (nbytes, sizeof (struct curv_points));
  257.       curv_head_in[in]->prev = NULL;
  258.       curv_tail_in[in] = curv_head_in[in];
  259.     }
  260.     else {
  261.       curv_tail_in[in]->next = (struct curv_points *) calloc (nbytes, sizeof (struct curv_points));
  262.       curv_tail_in[in]->next->prev = curv_tail_in[in];
  263.       curv_tail_in[in] = curv_tail_in[in]->next;
  264.     }
  265.     if (curv_head_out[out] == NULL) {
  266.       curv_head_out[out] = (struct curv_points *) calloc (nbytes, sizeof (struct curv_points));
  267.       curv_head_out[out]->prev = NULL;
  268.       curv_tail_out[out] = curv_head_out[out];
  269.     }
  270.     else {
  271.       curv_tail_out[out]->next = (struct curv_points *) calloc (nbytes, sizeof (struct curv_points));
  272.       curv_tail_out[out]->next->prev = curv_tail_out[out];
  273.       curv_tail_out[out] = curv_tail_out[out]->next;
  274.     }
  275.     ncp++;
  276.  
  277.     curv_tail_in[in]->same_point = curv_tail_out[out];
  278.     curv_tail_in[in]->next = NULL;
  279.     curv_tail_out[out]->next = NULL;
  280.     curv_tail_in[in]->xn = curv_tail_out[out]->xn = 2 * (j + 1);
  281.     curv_tail_in[in]->yn = curv_tail_out[out]->yn = 2 * (i + 1) + 1;
  282.     curv_tail_in[in]->edge_in = curv_tail_out[out]->edge_in = in;
  283.     curv_tail_in[in]->edge_out = curv_tail_out[out]->edge_out = out;
  284.   }
  285.   return (1);
  286. }
  287.  
  288.  
  289. /*
  290.  * link up all curvature points into closed polygons
  291.  */
  292. void
  293. linkage (Image * imgIn, Image * imgOut, int value)
  294. {
  295.   struct polygon *dummy_poly_ptr;
  296.   struct curv_points *first_point;
  297.   float *d_phi_head, *d_l_head;
  298.   int m;
  299.   int c;
  300.   char in_buf[IN_BUF_LEN];
  301.  
  302.   do {
  303.     d_phi_head = delta_phik;
  304.     d_l_head = delta_lk;
  305.     for (m = 0; m < NDIRECTIONS; m++)
  306.       if ((first_point = curv_head_out[m]) != NULL) {
  307.         poly (first_point);
  308.         break;
  309.       }
  310.   } while (first_point != NULL);
  311.  
  312.   if (ncp != 0)
  313.     printf ("...something wrong, didn't find all polygons!\n");
  314.   if (poly_head_ptr == NULL) {
  315.     printf ("...no objects found to analyze!\n");
  316.     return;
  317.   }
  318.  
  319. /*
  320.  * polygon analysis
  321.  */
  322.   printf ("...analyze a boundary?(y/n) -- ");
  323.   while ((c = readlin (in_buf)) == 'y') {
  324.     dummy_poly_ptr = select_poly (poly_head_ptr, imgOut, value);
  325.     printf ("\n...analyze another boundary?(y/n)");
  326.   }
  327. }
  328.  
  329. /*
  330.  * create a closed polygon given a starting curvature point 
  331.  */
  332. void
  333. poly (first_point)
  334.      struct curv_points *first_point;
  335. {
  336.   struct curv_points *next_point, *prev_point, *match_in_list ();
  337.  
  338.   if (poly_head_ptr == NULL) {
  339.     poly_head_ptr = (struct polygon *) calloc (nbytes, sizeof (struct polygon));
  340.     poly_tail_ptr = poly_head_ptr;
  341.     poly_head_ptr->prev_poly = NULL;
  342.   }
  343.   else {
  344.     poly_tail_ptr->next_poly = (struct polygon *) calloc (nbytes, sizeof (struct polygon));
  345.     poly_tail_ptr->next_poly->prev_poly = poly_tail_ptr;
  346.     poly_tail_ptr = poly_tail_ptr->next_poly;
  347.   }
  348.   poly_tail_ptr->next_poly = NULL;
  349.  
  350.   prev_point = first_point;
  351.   next_point = match_in_list (prev_point);
  352.   poly_tail_ptr->d_phi_ptr = delta_phik;
  353.   poly_tail_ptr->d_l_ptr = delta_lk;
  354.   poly_tail_ptr->poly_points = 0;
  355.   poly_tail_ptr->total_delta_phi = (float) 0;
  356.   poly_tail_ptr->first_x = first_point->xn;
  357.   poly_tail_ptr->first_y = first_point->yn;
  358.   poly_tail_ptr->first_edge_out = first_point->edge_out;
  359.   create_edge (prev_point, next_point);
  360.   delete_list (prev_point, OUT_LIST);
  361.   prev_point = next_point->same_point;
  362.   delete_list (next_point, IN_LIST);
  363.   ncp--;
  364.   poly_tail_ptr->poly_points++;
  365.   poly_tail_ptr->total_delta_phi = poly_tail_ptr->total_delta_phi + *delta_phik;
  366.   delta_lk++;
  367.   delta_phik++;
  368.   do {
  369.     next_point = match_in_list (prev_point);
  370.     create_edge (prev_point, next_point);
  371.     delete_list (prev_point, OUT_LIST);
  372.     prev_point = next_point->same_point;
  373.     delete_list (next_point, IN_LIST);
  374.     ncp--;
  375.     poly_tail_ptr->poly_points++;
  376.     poly_tail_ptr->total_delta_phi = poly_tail_ptr->total_delta_phi
  377.       + *delta_phik;
  378.     delta_lk++;
  379.     delta_phik++;
  380.  
  381.     if (ncp < 0) {
  382.       printf ("Something wrong, points<0\n");
  383.       exit (1);
  384.     }
  385.   }
  386.   while (prev_point != first_point);
  387.  
  388.   if (poly_tail_ptr->poly_points <= MIN_POLY_SIZE) {
  389.     if (poly_tail_ptr->prev_poly == NULL)
  390.       poly_head_ptr = NULL;
  391.     else {
  392.       poly_tail_ptr = poly_tail_ptr->prev_poly;
  393.       free (poly_tail_ptr->next_poly);
  394.       poly_tail_ptr->next_poly = NULL;
  395.     }
  396.   }
  397. }
  398.  
  399. /*
  400.  * given a point in the out-list, find the next point
  401.  * in the in-list with edge-out(out-list) = edge-in(in-list)
  402.  */
  403. struct curv_points *
  404. match_in_list (current_point)
  405.      struct curv_points *current_point;
  406. {
  407.   struct curv_points *tail_point, *head_point;
  408.   int direction;
  409.   int xcurr, ycurr, xnext, ynext;
  410.  
  411.   xcurr = current_point->xn;
  412.   ycurr = current_point->yn;
  413. #ifdef DPRINT
  414.   printf ("match_in_list: current point = (%d, %d)\n", xcurr, ycurr);
  415. #endif
  416.   direction = current_point->edge_out;
  417.   tail_point = curv_tail_in[direction];
  418.   head_point = curv_head_in[direction];
  419.   while ((head_point != NULL) || (tail_point != NULL)) {
  420.     switch (direction) {
  421.     case 0:
  422.       ynext = head_point->yn;
  423.       xnext = head_point->xn;
  424.       if ((ycurr == ynext) && (xnext > xcurr))
  425.         return (head_point);
  426.       break;
  427.     case 4:
  428.       ynext = tail_point->yn;
  429.       xnext = tail_point->xn;
  430.       if ((ycurr == ynext) && (xnext < xcurr))
  431.         return (tail_point);
  432.       break;
  433.     case 2:
  434.       ynext = tail_point->yn;
  435.       xnext = tail_point->xn;
  436.       if ((xcurr == xnext) && (ynext < ycurr))
  437.         return (tail_point);
  438.       break;
  439.     case 6:
  440.       ynext = head_point->yn;
  441.       xnext = head_point->xn;
  442.       if ((xcurr == xnext) && (ynext > ycurr))
  443.         return (head_point);
  444.       break;
  445.     case 1:
  446.       ynext = tail_point->yn;
  447.       xnext = tail_point->xn;
  448.       if (((xcurr + ycurr) == (xnext + ynext)) &&
  449.           (xnext > xcurr) && (ynext < ycurr))
  450.         return (tail_point);
  451.       break;
  452.     case 5:
  453.       ynext = head_point->yn;
  454.       xnext = head_point->xn;
  455.       if (((xcurr + ycurr) == (xnext + ynext)) &&
  456.           (xnext < xcurr) && (ynext > ycurr))
  457.         return (head_point);
  458.       break;
  459.     case 3:
  460.       ynext = tail_point->yn;
  461.       xnext = tail_point->xn;
  462.       if (((xcurr - ycurr) == (xnext - ynext)) &&
  463.           (xnext < xcurr) && (ynext < ycurr))
  464.         return (tail_point);
  465.       break;
  466.     case 7:
  467.       ynext = head_point->yn;
  468.       xnext = head_point->xn;
  469.       if (((xcurr - ycurr) == (xnext - ynext)) &&
  470.           (xnext > xcurr) && (ynext > ycurr))
  471.         return (head_point);
  472.       break;
  473.     }
  474.     tail_point = tail_point->prev;
  475.     head_point = head_point->next;
  476.   }
  477.   printf ("...no match found in outlist[%d]\n", direction);
  478.   exit (1);
  479. }
  480.  
  481. /*
  482.  * delete curvature point from either in-list or out-list
  483.  */
  484. void
  485. delete_list (current_point, list)
  486.      struct curv_points *current_point;
  487.      int list;
  488. {
  489.   int direction;
  490.  
  491.   if (current_point == NULL)
  492.     printf ("...current point can't be NULL!\n");
  493.  
  494.   if ((current_point->prev == NULL) && (current_point->next == NULL)) {
  495.     if (list == IN_LIST) {
  496.       direction = current_point->edge_in;
  497.       curv_head_in[direction] = NULL;
  498.       curv_tail_in[direction] = NULL;
  499.     }
  500.     else {
  501.       direction = current_point->edge_out;
  502.       curv_head_out[direction] = NULL;
  503.       curv_tail_out[direction] = NULL;
  504.     }
  505.   }
  506.   else if (current_point->prev == NULL) {
  507.     if (list == IN_LIST) {
  508.       direction = current_point->edge_in;
  509.       curv_head_in[direction] = current_point->next;
  510.     }
  511.     else {
  512.       direction = current_point->edge_out;
  513.       curv_head_out[direction] = current_point->next;
  514.     }
  515.     current_point->next->prev = NULL;
  516.   }
  517.   else if (current_point->next == NULL) {
  518.     if (list == IN_LIST) {
  519.       direction = current_point->edge_in;
  520.       curv_tail_in[direction] = current_point->prev;
  521.     }
  522.     else {
  523.       direction = current_point->edge_out;
  524.       curv_tail_out[direction] = current_point->prev;
  525.     }
  526.     current_point->prev->next = NULL;
  527.   }
  528.   else {
  529.     current_point->prev->next = current_point->next;
  530.     current_point->next->prev = current_point->prev;
  531.   }
  532.   free (current_point);
  533. }
  534.  
  535.  
  536. void
  537. print_curv ()
  538. {
  539.   int x, y;
  540.   struct curv_points *temp_head;
  541.  
  542.   for (x = 0; x < NDIRECTIONS; x++) {
  543.     y = 0;
  544.     temp_head = curv_head_in[x];
  545.     while (temp_head != NULL) {
  546.       y++;
  547.       temp_head = temp_head->next;
  548.     }
  549.     if (y != 0)
  550.       printf ("inlist[%d] = %d\n", x, y);
  551.   }
  552.  
  553.   for (x = 0; x < NDIRECTIONS; x++) {
  554.     y = 0;
  555.     temp_head = curv_head_out[x];
  556.     while (temp_head != NULL) {
  557.       y++;
  558.       temp_head = temp_head->next;
  559.     }
  560.     if (y != 0)
  561.       printf ("outlist[%d] = %d\n", x, y);
  562.   }
  563. }
  564.  
  565.  
  566. /*
  567.  * create an edge from old-point to new-point
  568.  */
  569. void
  570. create_edge (old_point, new_point)
  571.      struct curv_points *old_point, *new_point;
  572. {
  573.   int phi;
  574.  
  575.   switch (new_point->edge_in) {
  576.   case 0:
  577.   case 4:
  578.     *delta_lk = (float) abs (new_point->xn - old_point->xn);
  579.     phi = new_point->edge_out - new_point->edge_in;
  580.     if (phi == 7)
  581.       phi = -1;
  582.     *delta_phik = (float) phi;
  583.     break;
  584.   case 2:
  585.   case 6:
  586.     *delta_lk = (float) abs (new_point->yn - old_point->yn);
  587.     *delta_phik = (float) (new_point->edge_out - new_point->edge_in);
  588.     break;
  589.   default:
  590.     *delta_lk = (float) abs (new_point->xn - old_point->xn) * (float) ROOT2;
  591.     phi = new_point->edge_out - new_point->edge_in;
  592.     if (phi == -7)
  593.       phi = 1;
  594.     else if (phi == -6)
  595.       phi = 2;
  596.     else if (phi == 6)
  597.       phi = -2;
  598.     *delta_phik = (float) phi;
  599.     break;
  600.   }
  601. }
  602.